From 2e9aa6ca68b30facd3e4f29db78aa25d89a55c7e Mon Sep 17 00:00:00 2001 From: tsteven4 <13596209+tsteven4@users.noreply.github.com> Date: Sun, 9 Feb 2020 12:36:54 -0700 Subject: [PATCH] Rework interpolate filter, retire some magic numbers. (#498) * rewrite interpolate filter. What started out as an investigation of a clang-tidy warning about a floating point loop counter (FloatLoopCounter, cert-flp30-c) turned into this. Use integer loop counters. Use a variable step size, that is constant for each span, to meet the constraint of no interval > limit. Previously steps of the limit size were used leaving the last step in a span of a different likely shorter size. Be defensive about user input and track/route input, no more infinite loops for negative time limits. Support milliseconds. Reduce duplication. Enhance test to more easily verify interpolation steps. More efficient route list maniupulation. Correct documentation. Correct signed literals use with the unsigned argtype. This does reiterate, but not fix, the deficiency in route_head copying that recently came up where new members of route_head are not copied. * correct test to use long long format for TIMET_TIME_MS. * Eliminate some magic numbers for unit conversion. Previously the conversion factor for meters to feet was inaccurate. Correcting this caused a couple kml test cases to change. The ozi test was problematic, an altitiude was provided that endend in 0.5 feet, and was printed out rounded to feet. A one bit error in the mantissa of the double involved in the round trip to meters could change the rounding result. To fix this the ozi writer was changed to print track altitudes to a tenth of a foot. This matches an ozi track file found in the wild as well as our own ozi track sample. The classic-2 test involving gnuplot.style printed altitude to one millonth of a foot, i.e. 0.3 micrometers! The gnu plot style was changed to print to the thousandth of a foot, i.e. 0.3 millimeters. --- an1.cc | 4 +- arcdist.cc | 2 +- defs.h | 41 +- grtcirc.cc | 2 +- interpolate.cc | 197 +-- interpolate.h | 20 +- ozi.cc | 2 +- radius.cc | 2 +- reference/bounds-test.kml | 1500 ++++++++++---------- reference/dusky.gnuplot | 126 +- reference/gnuplot.style | 2 +- reference/track/20070813_short~plt.plt | 88 +- reference/track/bounds-test-track.kml | 1748 ++++++++++++------------ reference/track/interptrack.csv | 36 + reference/track/interptrack.gpx | 198 +-- reference/track/simpletrack.gpx | 2 +- reference/track/tinterptrack.csv | 12 + reference/track/tinterptrack.gpx | 20 +- saroute.cc | 2 +- testo.d/interpolate.test | 16 +- trackfilter.cc | 2 +- vecs.cc | 2 +- xmldoc/filters/interpolate.xml | 2 +- 23 files changed, 2061 insertions(+), 1965 deletions(-) create mode 100644 reference/track/interptrack.csv create mode 100644 reference/track/tinterptrack.csv diff --git a/an1.cc b/an1.cc index 7ecb0b8d2..af9c1d7f7 100644 --- a/an1.cc +++ b/an1.cc @@ -1167,11 +1167,11 @@ wr_init(const QString& fname) if (opt_zoom) { opt_zoom_num = atoi(opt_zoom); } - radius = .1609344; /* 1/10 mi */ + radius = .1609344; /* 1/10 mi in kilometers */ if (opt_radius) { radius = atof(opt_radius); if (!strchr(opt_radius,'k') && !strchr(opt_radius,'K')) { - radius *= 5280*12*2.54/100000; + radius *= kKilometersPerMile; } } } diff --git a/arcdist.cc b/arcdist.cc index edca0795c..5eee2a6ca 100644 --- a/arcdist.cc +++ b/arcdist.cc @@ -231,7 +231,7 @@ void ArcDistanceFilter::init() if ((*fm == 'k') || (*fm == 'K')) { /* distance is kilometers, convert to mile */ - pos_dist *= .6214; + pos_dist *= kMilesPerKilometer; } } } diff --git a/defs.h b/defs.h index ed732d202..b1a9fc486 100644 --- a/defs.h +++ b/defs.h @@ -66,18 +66,33 @@ # define M_PI 3.14159265358979323846 #endif -constexpr double FEET_TO_METERS(double feetsies) { return (feetsies) * 0.3048; } -constexpr double METERS_TO_FEET(double meetsies) { return (meetsies) * 3.2808399; } +/* + * The constants marked "exact in decimal notation" may be more accurately + * converted to doubles when we don't ask the compiler to do any multiplication. + * e.g. kMetersPerMile = 0.0254 * 12.0 * 5280.0; might have a larger error than + * kMetersPerMile = 1609.344; + * Historically we have had some problematic test cases where this mattered. + * It is a better idea to use test cases that aren't so demanding. + */ +constexpr double kMetersPerFoot = 0.3048; /* exact in decimal notation */ +constexpr double kFeetPerMeter = 1.0 / kMetersPerFoot; +constexpr double kMetersPerMile = 1609.344; /* exact in decimal notation */ +constexpr double kMilesPerMeter = 1.0 / kMetersPerMile; +constexpr double kKilometersPerMile = 1.609344; /* exact in decimal notation */ +constexpr double kMilesPerKilometer = 1.0 / kKilometersPerMile; + +constexpr double FEET_TO_METERS(double feetsies) { return (feetsies) * kMetersPerFoot; } +constexpr double METERS_TO_FEET(double meetsies) { return (meetsies) * kFeetPerMeter; } constexpr double NMILES_TO_METERS(double a) { return a * 1852.0;} /* nautical miles */ constexpr double METERS_TO_NMILES(double a) { return a / 1852.0;} -constexpr double MILES_TO_METERS(double a) { return (a) * 1609.344;} -constexpr double METERS_TO_MILES(double a) { return (a) / 1609.344;} +constexpr double MILES_TO_METERS(double a) { return (a) * kMetersPerMile;} +constexpr double METERS_TO_MILES(double a) { return (a) * kMilesPerMeter;} constexpr double FATHOMS_TO_METERS(double a) { return (a) * 1.8288;} -constexpr double CELSIUS_TO_FAHRENHEIT(double a) { return (((a) * 1.8) + 32);} -constexpr double FAHRENHEIT_TO_CELSIUS(double a) { return (((a) - 32) / 1.8);} +constexpr double CELSIUS_TO_FAHRENHEIT(double a) { return (((a) * 1.8) + 32.0);} +constexpr double FAHRENHEIT_TO_CELSIUS(double a) { return (((a) - 32.0) / 1.8);} constexpr long SECONDS_PER_HOUR = 60L * 60; constexpr long SECONDS_PER_DAY = 24L * 60 * 60; @@ -930,13 +945,13 @@ void setshort_whitespace_ok(short_handle, int n); void setshort_repeating_whitespace_ok(short_handle, int n); void setshort_defname(short_handle, const char* s); -#define ARGTYPE_UNKNOWN 0x00000000 -#define ARGTYPE_INT 0x00000001 -#define ARGTYPE_FLOAT 0x00000002 -#define ARGTYPE_STRING 0x00000003 -#define ARGTYPE_BOOL 0x00000004 -#define ARGTYPE_FILE 0x00000005 -#define ARGTYPE_OUTFILE 0x00000006 +#define ARGTYPE_UNKNOWN 0x00000000U +#define ARGTYPE_INT 0x00000001U +#define ARGTYPE_FLOAT 0x00000002U +#define ARGTYPE_STRING 0x00000003U +#define ARGTYPE_BOOL 0x00000004U +#define ARGTYPE_FILE 0x00000005U +#define ARGTYPE_OUTFILE 0x00000006U /* REQUIRED means that the option is required to be set. * See also BEGIN/END_REQ */ diff --git a/grtcirc.cc b/grtcirc.cc index 10846317a..d10e00369 100644 --- a/grtcirc.cc +++ b/grtcirc.cc @@ -57,7 +57,7 @@ static double dotproduct(double x1, double y1, double z1, double radtomiles(double rads) { - const double radmiles = EARTH_RAD * 100.0 / 2.54 / 12.0 / 5280.0; + const double radmiles = METERS_TO_MILES(EARTH_RAD); return (rads * radmiles); } diff --git a/interpolate.cc b/interpolate.cc index 9547688c3..ee4c7be41 100644 --- a/interpolate.cc +++ b/interpolate.cc @@ -19,15 +19,19 @@ */ -#include // for atoi, strtod +#include // for INT_MAX +#include // for abs, ceil, isfinite, round +#include // for abs, atoi, strtod #include // for QString -#include // for qAsConst, QAddConst<>::Type, foreach +#include // for qAsConst, QAddConst<>::Type #include "defs.h" #include "interpolate.h" #include "grtcirc.h" // for linepart, RAD, gcdist, radtomiles #include "src/core/datetime.h" // for DateTime +#include "src/core/logging.h" // for Fatal +#include "src/core/optional.h" // for optional #if FILTERS_ENABLED @@ -35,100 +39,111 @@ void InterpolateFilter::process() { - RouteList* backuproute = nullptr; - double lat1 = 0, lon1 = 0; - double altitude1 = unknown_alt; - unsigned int time1 = 0; - double frac; - - if (opt_route) { - route_backup(&backuproute); - route_flush_all_routes(); + RouteList backuproute; + if (opt_route != nullptr) { + route_swap(backuproute); } else { - track_backup(&backuproute); - route_flush_all_tracks(); + track_swap(backuproute); } - if (backuproute->empty()) { - fatal(MYNAME ": Found no routes or tracks to operate on.\n"); + if (backuproute.empty()) { + Fatal() << MYNAME ": Found no routes or tracks to operate on."; } - for (const auto* rte_old : qAsConst(*backuproute)) { - + for (const auto* rte_old : qAsConst(backuproute)) { + // FIXME: Allocating a new route_head and copying the members one at a + // time is not maintainable. When new members are added it is likely + // they will not be copied here! + // We want a deep copy of everything but with an empty WaypointList. route_head* rte_new = route_head_alloc(); rte_new->rte_name = rte_old->rte_name; rte_new->rte_desc = rte_old->rte_desc; - rte_new->fs = fs_chain_copy(rte_old->fs); + rte_new->rte_urls = rte_old->rte_urls; rte_new->rte_num = rte_old->rte_num; - if (opt_route) { + rte_new->fs = fs_chain_copy(rte_old->fs); + rte_new->line_color = rte_old->line_color; + rte_new->line_width = rte_old->line_width; + rte_new->session = rte_old->session; + if (opt_route != nullptr) { route_add_head(rte_new); } else { track_add_head(rte_new); } + + double lat1 = 0; + double lon1 = 0; + double altitude1 = unknown_alt; + gpsbabel::DateTime time1; bool first = true; - foreach (const Waypoint* wpt, rte_old->waypoint_list) { + for (const Waypoint* wpt : rte_old->waypoint_list) { if (first) { first = false; } else { - if (opt_interval && - wpt->creation_time.toTime_t() - time1 > interval) { - for (unsigned int timen = time1+interval; - timen < wpt->creation_time.toTime_t(); - timen += interval) { - auto* wpt_new = new Waypoint(*wpt); - wpt_new->SetCreationTime(timen); - wpt_new->shortname = QString(); - wpt_new->description = QString(); - - frac = (double)(timen - time1) / (double)(wpt->creation_time.toTime_t() - time1); - linepart(lat1, lon1, - wpt->latitude, wpt->longitude, - frac, - &wpt_new->latitude, - &wpt_new->longitude); - if (altitude1 != unknown_alt && wpt->altitude != unknown_alt) { - wpt_new->altitude = altitude1 + frac * (wpt->altitude - altitude1); - } - if (opt_route) { - route_add_wpt(rte_new, wpt_new); - } else { - track_add_wpt(rte_new, wpt_new); - } + gpsbabel_optional::optional timespan; + if (wpt->creation_time.isValid() && time1.isValid()) { + timespan = wpt->creation_time.toMSecsSinceEpoch() - + time1.toMSecsSinceEpoch(); + } + gpsbabel_optional::optional altspan; + if (altitude1 != unknown_alt && wpt->altitude != unknown_alt) { + altspan = wpt->altitude - altitude1; + } + + // How many points need to be inserted? + double npts = 0; + if (opt_time != nullptr) { + if (!timespan.has_value()) { + Fatal() << MYNAME ": points must have valid times to interpolate by time!"; } - } else if (opt_dist) { - double rt1 = RAD(lat1); - double rn1 = RAD(lon1); - double rt2 = RAD(wpt->latitude); - double rn2 = RAD(wpt->longitude); - double curdist = gcdist(rt1, rn1, rt2, rn2); - curdist = radtomiles(curdist); - if (curdist > dist) { - for (double distn = dist; - distn < curdist; - distn += dist) { - auto* wpt_new = new Waypoint(*wpt); - frac = distn / curdist; - wpt_new->SetCreationTime(frac * (wpt->creation_time.toTime_t() - time1) + time1); - wpt_new->shortname = QString(); - wpt_new->description = QString(); - linepart(lat1, lon1, - wpt->latitude, wpt->longitude, - frac, - &wpt_new->latitude, - &wpt_new->longitude); - if (altitude1 != unknown_alt && wpt->altitude != unknown_alt) { - wpt_new->altitude = altitude1 + frac * (wpt->altitude - altitude1); - } - if (opt_route) { - route_add_wpt(rte_new, wpt_new); - } else { - track_add_wpt(rte_new, wpt_new); - } - } + // interpolate even if time is running backwards. + npts = std::abs(timespan.value()) / max_time_step; + } else if (opt_dist != nullptr) { + double distspan = radtomiles(gcdist(RAD(lat1), + RAD(lon1), + RAD(wpt->latitude), + RAD(wpt->longitude))); + npts = distspan / max_dist_step; + } + if (!std::isfinite(npts) || (npts >= INT_MAX)) { + Fatal() << MYNAME ": interpolation interval too small!"; + } + + // Insert the required points + int nmax = static_cast(ceil(npts)) - 1; // # of points to insert + for (int n = 0; n < nmax; ++n) { + double frac = static_cast(n + 1) / + static_cast(nmax + 1); + // We create the inserted point from the Waypoint at the end of the + // span. Another choice would be the Waypoint at the beginning of + // the span. We clear some fields but use a copy of the rest or the + // interpolated value. + auto* wpt_new = new Waypoint(*wpt); + wpt_new->shortname = QString(); + wpt_new->description = QString(); + if (timespan.has_value()) { + wpt_new->SetCreationTime(0, time1.toMSecsSinceEpoch() + + round(frac * timespan.value())); + } else { + wpt_new->creation_time = gpsbabel::DateTime(); + } + linepart(lat1, lon1, + wpt->latitude, wpt->longitude, + frac, + &wpt_new->latitude, + &wpt_new->longitude); + if (altspan.has_value()) { + wpt_new->altitude = altitude1 + (frac * altspan.value()); + } else { + wpt_new->altitude = unknown_alt; + } + if (opt_route != nullptr) { + route_add_wpt(rte_new, wpt_new); + } else { + track_add_wpt(rte_new, wpt_new); } } } - if (opt_route) { + if (opt_route != nullptr) { route_add_wpt(rte_new, new Waypoint(*wpt)); } else { track_add_wpt(rte_new, new Waypoint(*wpt)); @@ -137,31 +152,35 @@ void InterpolateFilter::process() lat1 = wpt->latitude; lon1 = wpt->longitude; altitude1 = wpt->altitude; - time1 = wpt->creation_time.toTime_t(); + time1 = wpt->creation_time; } } - backuproute->flush(); - delete backuproute; + backuproute.flush(); } void InterpolateFilter::init() { - char* fm; - if (opt_interval && opt_dist) { - fatal(MYNAME ": Can't interpolate on both time and distance.\n"); - } else if (opt_interval && opt_route) { - fatal(MYNAME ": Can't interpolate routes on time.\n"); - } else if (opt_interval) { - interval = atoi(opt_interval); - } else if (opt_dist) { - dist = strtod(opt_dist, &fm); + if ((opt_time != nullptr) && (opt_dist != nullptr)) { + Fatal() << MYNAME ": Can't interpolate on both time and distance."; + } else if ((opt_time != nullptr) && (opt_route != nullptr)) { + Fatal() << MYNAME ": Can't interpolate routes on time."; + } else if (opt_time != nullptr) { + max_time_step = 1000 * strtod(opt_time, nullptr); // milliseconds + if (max_time_step <= 0) { + Fatal() << MYNAME ": interpolation time should be positve!"; + } + } else if (opt_dist != nullptr) { + max_dist_step = strtod(opt_dist, &fm); if ((*fm == 'k') || (*fm == 'K')) { /* distance is kilometers, convert to miles */ - dist *= .6214; + max_dist_step *= kMilesPerKilometer; + } + if (max_dist_step <= 0) { + Fatal() << MYNAME ": interpolation distance should be positve!"; } } else { - fatal(MYNAME ": No interval specified.\n"); + Fatal() << MYNAME ": No interval specified."; } } diff --git a/interpolate.h b/interpolate.h index 36fa127de..c3ba39add 100644 --- a/interpolate.h +++ b/interpolate.h @@ -22,10 +22,12 @@ #ifndef INTERPOLATE_H_INCLUDED_ #define INTERPOLATE_H_INCLUDED_ -#include // for QVector +#include // for QVector +#include // for qint64 -#include "defs.h" // for ARG_NOMINMAX, arglist_t, ARGTYPE_BEGIN_EXCL, ARG... -#include "filter.h" // for Filter +#include "defs.h" // for ARG_NOMINMAX, arglist_t, ARGTYPE_BEGIN_EXCL, ARG... +#include "filter.h" // for Filter +#include "src/core/optional.h" // for optional #if FILTERS_ENABLED @@ -40,15 +42,15 @@ public: void process() override; private: - char* opt_interval = nullptr; - unsigned int interval = 0; - char* opt_dist = nullptr; - double dist = 0; - char* opt_route = nullptr; + char* opt_time{nullptr}; + double max_time_step{0}; + char* opt_dist{nullptr}; + double max_dist_step{0}; + char* opt_route{nullptr}; QVector args = { { - "time", &opt_interval, "Time interval in seconds", nullptr, + "time", &opt_time, "Time interval in seconds", nullptr, ARGTYPE_BEGIN_EXCL | ARGTYPE_BEGIN_REQ | ARGTYPE_INT, "0", nullptr, nullptr }, diff --git a/ozi.cc b/ozi.cc index 184b86acd..e309f3cda 100644 --- a/ozi.cc +++ b/ozi.cc @@ -309,7 +309,7 @@ ozi_track_disp(const Waypoint* waypointp) *stream << qSetRealNumberPrecision(6) << waypointp->latitude << ',' << waypointp->longitude << ',' << new_track << ',' - << qSetRealNumberPrecision(0) << alt << ',' + << qSetRealNumberPrecision(1) << alt << ',' << ozi_time << ",,\r\n"; new_track = 0; diff --git a/radius.cc b/radius.cc index 37588fcf2..c52389911 100644 --- a/radius.cc +++ b/radius.cc @@ -145,7 +145,7 @@ void RadiusFilter::init() if ((*fm == 'k') || (*fm == 'K')) { /* distance is kilometers, convert to miles */ - pos_dist *= .6214; + pos_dist *= kMilesPerKilometer; } } diff --git a/reference/bounds-test.kml b/reference/bounds-test.kml index 1cae4cdff..65298dd15 100644 --- a/reference/bounds-test.kml +++ b/reference/bounds-test.kml @@ -7,7 +7,7 @@ 36.438270 87257.203848 - + - + - + - + - +